function [ Fxs,Fys ] = Bearingforce(ts,w,xs,ys,vxs,vys)
%  x,y: displacements of inner ring in X and Y directions
%  vx,vy: velocities of inner ring in X and Y directions
%  Kb: contact deformation coefficient between roller and raceway
%  di,de: inner ring diameter, outer ring diameter
%  w: ratating speed of inner ring, it is set to zero here.
%  output: rolling bearing forces, Fx, Fy [matrix version]
	Kb=1.8e8;	
	Kn=1.0831e10;
	Cb=10e-5*Kb;                                                                % damping coefficient
	ri=0.023;              
	re=0.031;              
	Nb=16;                                                                      %roller number
	wc=ri*w/(ri+re);                                                            %cage speed
	gama0=10e-6;                                                                %radial clearance
	Fxs=zeros(size(xs));
	Fys=zeros(size(ys));                    

	for i=1:Nb
		theta=2*pi*(i-1)/Nb+wc*ts;                                              % column vector
		tempvalue=xs.*sin(theta)+ys.*cos(theta)-gama0;                          % matrix
		tempvalue(tempvalue<0)=0;
		Fxs=Fxs-Kn*(tempvalue.^(3/2)).*sin(theta);
		Fys=Fys-Kn*(tempvalue.^(3/2)).*cos(theta);   
	end
	Fxs=Fxs-Cb*vxs;                                                             % considering linear damping in x direction
	Fys=Fys-Cb*vys;                                                             % considering linear damping in y direction
end